#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;

/****************************************************************
 * Homework assignment 2, CSCE 240, Spring 2010
 * Author: Duncan Buell
 * Date last modified: 10 January 2010
 *
 * This program computes an approximation for 'pi' using the
 * standard Monte Carlo approach.
 *
 * The program reads from standard input the number of iterations
 * to perform and a "print frequency".  It then calls a function
 * that runs a loop for that many iterations, generating random
 * numbers 'x' and 'y' to be taken as '(x,y)' coordinates inside
 * a unit square in the northeast quadrant from '(0,0)'.  If the
 * point 'P = (x,y)' lies inside the unit circle, we count it
 * toward the value of 'pi'.  If it lies outside the circle, it
 * does not count.
 *
 * To monitor the progress of the function, we print out a
 * formatted line every 'printFrequency' iterations.
 *
 * NOTE:  In a real computation, we would not run a specific
 * number of iterations but would instead run a 'while' loop
 * until the differences between the 'old' approximation and
 * a 'new' approximation fell below some small value.
**/

double computePi(int numberOfIterations, int printFrequency);

int main( )
{
  int numberOfIterations; // The number of iterations to perform.
  int printFrequency;     // How often to print the current approx.
  double approxPi;        // The approximation computed.

  cout << "Enter the number of iterations: "; 
  cin >> numberOfIterations;
  cout << "Enter the frequency with which to print the current approx: "; 
  cin >> printFrequency;
  cout << "We approximate 'pi' with " << numberOfIterations
       << " iterations." << endl;
  cout << "We print the current approximation every "
       << printFrequency
       << " iterations." << endl;

  if(numberOfIterations <= 0)
  {
    cout << "ERROR: illegal nonpositive value "
         << numberOfIterations
         << " for 'numberOfIterations'.\n";
    exit(0);
  }
  if(printFrequency <= 0)
  {
    cout << "ERROR: illegal nonpositive value "
         << printFrequency
         << " for 'printFrequency'.\n";
    exit(0);
  }

  approxPi = computePi(numberOfIterations, printFrequency);

  cout << "The approximate value of 'pi' is " << approxPi << endl;

  return 0;
} // int main( )

/****************************************************************
 * Function to compute an approximation for 'pi' using the
 *   standard Monte Carlo method.
 * Author: Duncan Buell
 * Date last modified: 2 January 2010
 *
 * This program computes an approximation for 'pi' using the
 * standard Monte Carlo approach.
 *
 * The program takes as arguments the number of iterations
 * to perform and a frequency with which to print intermediate
 * approximations.  It then runs a loop for that many iterations,
 * generating random numbers 'x' and 'y' to be taken as '(x,y)'
 * coordinates inside a unit square in the northeast quadrant
 * from '(0,0)'.  If the point 'P = (x,y)' lies inside the unit
 * circle, we count it toward the value of 'pi'.  If it lies
 * outside the circle, it does not count.
 *
 * A point is inside the circle if the length of the line from
 * '(0,0)' to '(x,y)' is less than or equal to 1.0, and outside
 * otherwise. 
 *
 * Input:
 *   numberOfIterations: the 'int' number of iterations in
 *                       the loop we will execute
 *   printFrequency: the frequency with which to print the
 *                   current approximation
 * Output:
 *   This function prints to 'cout' a formatted table every
 *   'printFrequency' iterations.
 *
 * NOTE: We print a line when 'iteration' modulo 'printFrequency'
 *       is zero.  We also print an initial zero value for the
 *       zeroth iteration just for completeness.
 *
 * NOTE: We assume that the two input arguments are positive.
**/
double computePi(int numberOfIterations, int printFrequency)
{
  // Constants in lieu of magic numbers.
  // Note that we need the const variable for 'RANDMAXDOUBLE'
  // since 'enum' items are only of 'int' type.
  enum {ITERATIONPRINTWIDTH = 8}; // field width for output
  enum {PIPRINTPRECISION = 8};    // # of decimals for output
  const double RANDMAXDOUBLE = static_cast<double>(RAND_MAX);

  int iteration; // Iteration number.
  int count; // The number of times we have generated a point
             // inside the unit circle.
  double pi;   // Our computed approximation to 'pi'
  double x, y; // The 'x' and 'y' coordinates.
  double r;    // The length of the vector from '(0,0)' to '(x,y)'.

  // We initialize values here so we can have our iteration
  // zero output line match up in format with the rest from
  // inside the loop.  This is really the only reason we didn't
  // declare the 'iteration' variable in the 'for' statement. 
  iteration = 0;
  count = 0;
  pi = 0.0;

  // And then output for iteration zero before the loop begins.
  cout.setf(ios::right);
  cout.width(ITERATIONPRINTWIDTH);
  cout << iteration << " ";
  cout.setf(ios::right);
  cout.width(ITERATIONPRINTWIDTH);
  cout << count << " ";
  cout.setf(ios::fixed);
  cout.setf(ios::showpoint);
  cout.precision(PIPRINTPRECISION);
  cout << pi << endl;

  // Now run the loop.
  for(iteration = 1; iteration <= numberOfIterations; ++iteration)
  {
    // Generate random integers and divide down to force them into
    // the range 0.0 to 1.0.
    x = rand() / RANDMAXDOUBLE;
    y = rand() / RANDMAXDOUBLE;

    // Compute the length of the vector from '(0,0)' to '(x,y)'
    // and score this iteration as inside the unit circle if it
    // is in fact inside.
    r = sqrt(x*x + y*y);
    if(r <= 1.0)
    {
      ++count;
      pi = 4.0 * static_cast<double>(count) /
                   static_cast<double>(iteration);
    }

    // Output status every 'printFrequency' iterations.
    if(0 == iteration % printFrequency)
    {
      cout.setf(ios::right);
      cout.width(ITERATIONPRINTWIDTH);
      cout << iteration << " ";
      cout.setf(ios::right);
      cout.width(ITERATIONPRINTWIDTH);
      cout << count << " ";
      cout.setf(ios::fixed);
      cout.setf(ios::showpoint);
      cout.precision(PIPRINTPRECISION);
      cout << pi << endl;
    }
  } // for(int iteration = 1;...

  return pi;
} // double computePi(int numberOfIterations, int printFrequency)

